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Abstract: We non-perturbatively determine the renormalization constant and the im- 
provement coefficients relating the renormalized current and subtracted quark mass of 
(quenched) valence quarks propagating in a sea of 0(a) improved two massless quarks. 
We employ the Schrodinger functional scheme and fix the physical extent of the box by 
working at a constant value of the renormalized coupling. Our calculation yields results 
which cover two regions of bare parameter space. One is the weak-coupling region suitable 
for volumes of about half a fermi. By making simulations in this region, quarks as heavy 
as the bottom can be propagated with the full relativistic QCD action and renormalization 
problems in HQET can be solved non-perturbatively by a matching to QCD in finite vol- 
ume. The other region refers to the common parameter range in large-volume simulations 
of two-flavour lattice QCD, where our results have particular relevance for charm physics 
applications. 
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1. Introduction 

Hadronic observables can be calculated with non-perturbative accuracy by recurring to lat- 
tice QCD. Among the shortcomings of all the lattice formulations of the continuum theory 
are cutoff effects, an unavoidable source of systematic errors that is particularly relevant 
in the case of flavoured quantities. On the one hand, cutoff effects can be quantified and 
removed by repeating the calculations with different lattice spacings and by extrapolating 
the results to the continuum limit. On the other hand, the bigger are the leading cutoff 
effects the harder will be the extrapolations and/or the computational effort required to 
perform simulations at smaller values of the lattice spacing. 

There exist lattice formulations of the continuum theory carrying leading cutoff effects 
proportional to the square of the lattice spacing. 0(a) improvement is achieved at the 
price of a high computational cost by working with lattice Dirac operators, solutions of 
the celebrated Ginsparg- Wilson relation, that preserve at fixed cutoff a lattice version of 
chiral symmetry (see ref. [1] for a recent review). Automatic 0(a) improvement of some 
or all physical observables can be also achieved by introducing in the naive lattice fermion 
action a dimension-five operator preserving a non-singlet chiral symmetry while breaking 
at the same time the doubling symmetry and another symmetry of the continuum action, 
e.g. parity in the case of the twisted mass formulation [2]. 

In this paper we concentrate on the Wilson formulation of lattice QCD that is very 
well understood, provides a strictly local fermion action, preserves parity and all the vector 
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symmetries of the continuum theory and where 0(a) improvement is not automatic, but can 
be obtained by implementing Symanzik's program [3, 4, 5]. Furthermore, improvement can 
be restricted to on-shell quantities and correlation functions at physical distances without 
any loss [6]. 

In particular we calculate non-perturbatively with two massless dynamical flavours of 
sea quarks the coefficients 6m and 6a — 6p and the renormalization constant Z that are 
needed in order to renormalize (quenched) valence quark masses in the improved theory. 
The same quantities have been obtained in the quenched case in refs. [7, 8, 9], and here we 
closely follow the strategy devised in these papers. We employ the Schrodinger functional 
scheme [10, 11, 12, 13] and fix the physical extent of the box by working at a constant 
value of the renormalized coupling. 

Our results span two regions of bare parameter space. The first is a weak coupling 
region corresponding to a fixed physical volume of about 0.5 fm. Simulations on such a 
small volume are required in order to study heavy-light hadronic systems on the lattice with 
quarks as heavy as the bottom. Indeed, given the computing resources available nowadays, 
b-quarks can be propagated with their full relativistic action and without introducing 
unacceptably high cutoff effects only in sufficiently small volumes. These simulations can 
then be used to match non-perturbatively HQET with QCD [14] , and physical results can be 
obtained by calculating in the effective theory the running of the renormalization constants 
down to scales corresponding to sufficiently large physical volumes. Alternatively, finite- 
volume QCD results can be corrected by calculating in the full theory the finite-volume 
effects by means of a finite-size scaling technique [15] . These two strategies can be profitably 
combined as it has been done in the quenched theory [16]. 

The other region of parameter space of our results spans the common range of bare cou- 
plings used in large- volume simulations of two-flavour lattice QCD, where non-perturbative 
improvement has particular relevance for charm physics applications [17, 18, 19]. 

2. Improvement conditions at constant physics 

The basic idea of on-shell improvement of a lattice field theory such as Wilson's formulation 
of lattice QCD as employed in the present paper consists in adding counterterms to the 
unimproved lattice action and fields such that the leading cutoff effects are cancelled in 
on-shell amplitudes. Having identified the cutoff effect at fixed 50 1 where a continuum 
symmetry is broken, the coefficients of the different improvement terms need to be fixed by 
suitable improvement conditions. These conditions can either require pure lattice artifacts, 
i.e. combinations of observables that are known to approach zero in the continuum limit 
of the theory, to vanish, or even more generally the scaling of some renormalized quantity. 
This then defines the values of the improvement coefficients as a function of the lattice 
spacing a or, equivalently, of the bare gauge coupling, = 6//3. 

In perturbation theory, where in principle lattice artifacts can be obtained from any 
renormalized quantity by subtracting its value in the continuum limit, the improvement 
coefficients are unique. Beyond perturbation theory, where one wants to determine the 
improvement coefficients by Monte Carlo calculations, they do not remain unique and 
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depend on the choice of the condition. This in turn can be exploited in order to choose an 
improvement condition with a good sensitivity to the coefficient in question. In Wilson's 
lattice QCD, the resulting ambiguity is of 0(a) in the coefficients and of O(a^) in the 
physical observables after improvement. Extra care has to be exercised, however, if these 
ambiguities are large. A typical example of such a case, which has already been identified 
in the quenched approximation [8, 9], is the renormalization of the subtracted bare quark 
mass in the 0(a) improved theory, which amounts to consider QCD with non-degenerate 
quark masses. 

The authors of ref. [21] have discussed the details of improvement and renormalization 
of the quark field bilinears in the presence of mass non-degenerate quarks. According to 
their results a fiavour non-singlet bilinear O*-' renormalizes as 



O'^ = Zo 



l + 6o2(«"iq,, +amqj) + 6oiVfa"if*'''^ O'^ , (2.1 



while, because of the presence of disconnected diagrams, flavour singlet bilinears require 
a different renormalization constant and additional improvement coefficients. In the rest 
of our paper we consider only fiavour non-singlet quantities and work at vanishing sea 
quark masses so that we need neither to write down contributions proportional to the 
improvement coefficient bp nor to discuss renormalization of fiavour singlet quantities any 
further. 

Let us illustrate this now in more detail by directly entering the discussion of the 
situation in two-fiavour QCD addressed by this work. To avoid notational overload and 
since any renormalization scale dependence is unimportant for the purpose of this section, 
we write as shorthand for rn{fj,) in the following. 

2.1 The PCAC relation and mass non-degenerate quarks 

For the formulation of proper improvement conditions in a non-perturbative calculation 
it is generally advantageous to use lattice artifacts that derive from a symmetry of the 
continuum field theory that is not respected by the lattice regularization. A particularly 
suited choice for the purpose of this work are violations of the PCAC relation between the 
axial vector current and the pseudo-scalar density. Assuming that on-shell 0(a) improve- 
ment has been fully implemented and using renormalized currents in the lattice regularized 
theory, it reads 

( d^iA^r^ix) O ) = 2mn{ (Pnnx) O) , d^ = l{d, + d;) , (2.2) 

up to O(a^) corrections, as long as x is not in the support of the source operator O. 
We adopt the standard renormalization pattern of improved lattice currents in a mass 
independent scheme which is conventionally decomposed as 

_ ZA(l + 6Aamq) , 1 {d^{Ai)-ix) O) 

"^^-Zp(l + 5pamq)"^ + °(")' "^-2 (P«(x)0) " ^^'^^ 

Here, m denotes the bare current quark mass defined through some matrix element of the 
bare PCAC relation, where (^i)JI^ = + ac^d^P'^ is the 0(a) improved axial current and 
CA(ffo) is taken to be the non-perturbative value of the two-fiavour theory [20]. 
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On the other hand, starting from the bare quark mass mg of the Lagrangian, we have 
an alternative representation of the renormalized quark mass in terms of the 0(a) improved 
subtracted bare quark mass, fhq, viz. 

mR = Zm m-q , mq = mq(l + bmarriq) , mq = mo - m-c . (2.4) 

By equating eqs. (2.3) and (2.4) one obtains the bare current quark mass in terms of the 
subtracted bare quark mass: 

m = ^ (l + ^mamq)(l + 6pamq) ^ 
Za (1 + ftAomq) 



ZmZp 



Za 

The renormalization factor 



1 + (&m + &P - ^a) amq + 0([amq]2)j mq + 0{a'^) . (2.6) 



y(~2-. _ ^m(go,a/^)Zp(gg,a/i) 



is finite, because the divergence of Zm is canceled by that of Z-p. Therefore, Z is only a 
function of the improved bare coupling, 'qq = (7g(l + bgam!^^^^), while the second factor 
in (2.5) is not. Our simulations of the two-flavour theory were performed at vanishing sea 
quark mass (amq ~ 0) and the valence quarks entering the bilinears are treated in the 
quenched approximation, hence = holds throughout.^ To keep notation short, we 
do not explicitly write out the ^Q-dependence in our improvement coefficients and Z from 
now on. 

To disentangle Z, bm and 6p — 6a in eq. (2.6) from each other later on, one has to con- 
sider non-degenerate current quark masses, mR^j ^ '^Rj- t^^is theory the renormalized, 
improved axial currents and densities are given by a straightforward generalization of what 
was written down earlier [7]. More precisely, since the presence of non-degenerate masses 
breaks flavour symmetry, our starting point now are the off-diagonal bilinear fields, 

4^'(x) = V^,(x)7m75V',(x), P'\x) = ^,{x)^^i^j{x). (2.8) 

As above, the renormalized PCAC relation 

( d^{A^Y^{x) O ) = (mR,i + mR,,) ( (Pr)^^(x) O ) + 0{a^) (2.9) 

has to be understood as an operator identity, with 

{A^Y\x) = ZA{gl) [ 1 + 6a(<7o) 2 (a"iq,i + am^^j) ] , (2.10) 

(Pr)'^(x) = Zp(<72, a/i) [ 1 + 6p(5o) \ (amq,i + amq,,) ] P'\x) (2.11) 

and {Aiffl = A^fi + acAd^P^^ . The renormalization constants Za and Zp involved are the 
same as those in the theory with two mass-degenerate quarks, and the 6-coefficients now 



^For a more general treatment without these simphfications, the reader may consuh ref. [21]. 
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multiply cutoff effects proportional to the average ^ {aniq^i + arriq^j ) of the subtracted bare 
quark masses. 

Correspondingly, also the subtracted quark masses of each flavour i get separately 
improved, but share the renormalization factor of the degenerate case, viz. 

rriR^i = Zm rhq^i . fhq^i = ruq^i (1 + ftmaw-q.i) , ruq^i = mo^i - rric , (2.12) 

and in analogy to eq. (2.5), 

"iR,j = ^ Zniq^i (1 + bmarriq^i) , Z = "^'"^^ . (2.13) 

Zp Za 

By inserting the last equation into the PCAC relation (2.9), and by solving for the associ- 
ated bare current quark masses, we arrive at the generalization of eqs. (2.5) and (2.6) for 
the case of mass non-degenerate quarks: 

<(6»o4 + acAaaP*^)_0> ^ 



2 {P'i O) 



Z 



I (mq,i + m-qj) + I bra (om^ ^ + amlj) - j (6a - bp) a {m^^i + m^^jf 



again up to terms of O(a^). Moreover, we introduce a third renormalized mass defined as 
the mean value of two given ones, 

^ ("T-R,! +™'R,2) = -^m ("lq,l -hmq,2) + f'ma^ ('IT'q,! +"1^,2)] i (2-15) 

which differs from the one obtained by the renormalization of the mean mo, 3 = 2(^^0,1 + 
mo,2) taken at the level of bare quark masses, 

"IR,3 = Zra ?7lq,3(l + braam^^'^) = Z^ [I (^q^i + mq,2) + V O i (w-q,! + "lq,2)^ ] • (2.16) 

The difference of the quadratic terms in eq. (2.15) and (2.16) can be used to estimate 6m. 
By building appropriate differences with non-degenerate quarks it is possible to construct 
estimators for the improvement coefficient b^ and 6a — 6p separately. 

In the following we show how one can extract 6a — 6p, b^ and Z by applying eq. (2.14) 
with a few external operators resp. kinematical conditions [7, 8] (see also [22]). 



2.2 Schrodinger functional correlation functions 

Particularly in the context of the QCD Schrodinger functional (SF) and its lattice for- 
mulation [10, 11, 12, 13], the use of PCAC has proven to be very advantageous to solve 
improvement and renormalization problems in many situations, because there one has a 
large flexibility to choose appropriate improvement and renormalization conditions. Be- 
yond that, one can compute, e.g., the improvement coefficients also for small values of the 
bare coupling go, thereby making contact with their perturbative expansions. An overview 
is given, for instance, in ref. [23]. 
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If C and C denote boundary quark and anti-quark fields at Euclidean time xq = [24], 
the field products summed over the spatial volume, 

^'' = Z^ECH75 2^X(v) , 0"^ = ^J2^'{u)^,^^r^C'i^) , (2.17) 

u,v u,v 

yield pseudo-scalar boundary sources projected onto zero momentum. From these, the 
boundary-to-bulk resp. forward SF correlation functions in the pseudo-scalar channel are 
constructed from the axial current and density as 

/a(xo) = -y E (^o(^) O") , /p(^o) = -f E (^"(^) ^'') ■ (2.18) 

X X 

The corresponding bulk-to-boundary resp. backward correlation functions, gA,p, are related 
to the former by a time reflection and fall off with T — xq from their sources at xq = T as do 
the forward correlation functions with xq from xq = 0. In a vanishing background field, this 
time reflection symmetry can be utilized by taking their mean before building secondary 
quantities, in order to reduce the statistical noise. In the theory with mass non-degenerate 
quarks the bare correlation functions generalize to 

flM = -f E i^oi^) > , f^\xo) = -y E (^^'(^) > • (2.19) 

X X 

We do not apply a summation convention here, and the boundary sources O^^ are formally 
composed by non-degenerate boundary quark and anti-quark fields in the same way as the 
currents in eq. (2.8). Backward correlation functions g^^p{T — xq) are defined analogously. 

Including the 0(a) correction of the axial current as in eq. (2.10), the improved bare 
current quark masses^ expressed in terms of the just introduced SF correlation functions 
are given by 

/ x_ / 9o/a(3^o) + acA<9o5o/p (xo) , . 

mij{xo) = mij{xo;L/a,T/L,9) = ^ -j- — ^ ^ , (2.20) 

2/p'(a;o) 

where, only here, we explicitly indicate their additional dependence on L/a, T/L and the 
periodicity angle 9 in the boundary conditions of the fermion fields. In the degenerate case 
(i = j), niij reduces to the non-singlet PCAC mass of a fiavour degenerate doublet. In 
this definition, Oq denotes the average of the usual forward and backward derivatives but, 
as in [7, 8, 9], we have also employed current quark masses involving improved derivatives, 
obtained by the replacements 

do^do{l-la^dtdo) , dtdo ^ d^do {l - ^ a^dtdo) , (2.21) 

in order to eventually estimate &m, &a — bp and Z. When acting on smooth functions, these 
derivatives have an error of ©((^qO^, a^). 



^This expression for the PCAC masses is only 0(a) improved up to a factor 1 + 1 (6a — 6p) (arriq^i +c 
for quark mass dependent cutoff effects. 
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2.3 Estimators for 6a — bp, and Z 

Following the strategy already proposed for the quenched computation in [8, 9], we form 
ratios of suitable combinations of degenerate and non-degenerate current quark masses. 
Representing eq. (2.14) through the niij now, it turns into 



Z 



4 im„_i + mqj) + I 6m (a?«q,i + "?"qj) - i (&A - bp) a {niq^i + mqj)^ . (2.22 



This enables us to derive direct estimators for 6a — bp, 6ni and Z in terms of various bare 
PC AC quark masses ruij, which can by computed via SF correlation functions according 
to eq. (2.20). By exploiting the symmetry properties [8] of eq. (2.22) and keeping only the 
relevant terms, i.e. 

2ami2 - {aniii + 077122) = Z ^ {bA - bp) (amq^i - mq^2)^ [ 1 + 0(amq_i-|-amq,2)] , (2.23) 
amii - am22 = -Z' (amq,i — amq,2) [ 1 + 0(amq,i-|-amq,2)] , (2.24) 
ami2 - am^-i = \bmZ (amq^i - amq_2)^ [ 1 + 0(amq_i-|-amq^2) ] , (2.25) 

it is straightforward to deduce the following ratios, which define estimators for the desired 
improvement coefficients and renormalization constant: 

_ 2 (2mi2 - mil - 77122) u \ , , M /o or^ 

-Rap = 7 TT r = (6a - 6p) 1 + OfarTiq.i +a77iq,2) , (2.26) 

(77711 - 77122) (a"lq,l - amq,2) L v . 

R^=- 4 (777 1 2 - ^33) 6m[l + 0(a777q,i+a777q,2)] , (2.27) 

(77711 - 77722) (a77lq,l - a7n,q,2) 

Rz = ""^^ ~ "^^^ + {Rap - Rm) (amn + 077722) = Z + 0(a^) . (2.28) 

"7q,l - mq,2 

While the leading cutoff effects on the estimators for the improvement coefficients Rap and 
i?m depend (linearly) on the bare subtracted quark masses, the leading correction to the 
estimator Rz for the renormalization factor is O(a^) after inserting 

_ 2 (277133 - 77711 - 777 22) fi. u u \ , r^/ , \ fo on\ 

RAP-Rm = 7 TT r = (6A-0p-0m) + 0(a7r7q,i + a777q,2) . (2.29) 

(77711 - ^22) (a777q,l - a777q,2) 

The difference i?AP — Rm does not depend on the off-diagonal mass 77712 and the cutoff- 
dependence of Rz is the one expected for a renormalization constant in the 0(a) improved 
theory. Note that the essential piece for the disentanglement of -Rap and Rm is the intro- 
duction of a third mass in accordance with 

mo,3 = 5 {mo,l + 7770,2) ^ ?T7q,3 = 5 ("7q,l -|- 777q,2) , (2.30) 

as already anticipated in eqs. (2.15) and (2.16) before. 

Since the bare current quark masses 777 jj are functions of the time coordinate xq, 
where the operator insertion of the currents ^0 and P takes place, so do the estimators 
Rx, X = AP,m, Z, depend on xq. Therefore, one has to specify the exact timeslice 
range used to compute the estimators Rx- This choice must be considered as being part 
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of a particular improvement condition, which uniquely defines the ratios Rx as smooth 
functions of qq. Within the SF formalism, a good plateau in the local PCAC masses usually 
sets in around the midpoint of the time extension, T/2. This need not necessarily be the 
case for observables derived from these masses, especially when investigating improvement 
coefficients which genuinely emerge only at order a. Fluctuations are in general expected 
to show up, depending on the quantity and the region in bare parameter space where the 
simulation takes place, but this does not invalidate the once chosen improvement condition. 

As we will use T = 3L/2 later on, our final estimators are defined as timeslice averages 
over the central third of the time extension: 

(T-L/2)/a 

Rx^ ^ " ^ E Rx{xo) , XG{AP,m,Z}. (2.31) 

1 — LI + a ^-^ 

' xo=L/2a 

Note that the main reason for this is rather to be less dependent on a special timeslice 
— which means to get a smoother and more reliable signal — than simply decreasing 
the size of the error. Furthermore, since the number of timeslices used for the plateau 
average scales with the lattice size, the physical size of this plateau is kept constant and 
hence automatically respects our requirement of a constant physics condition to be used 
below. subsectionGeneral remarks on improvement conditions at constant physics Before 
proceeding with the actual computation of 6a — 6p, and Z, we would like to discuss the 
virtues of improvement conditions at constant physics on more general grounds [23]. 

In the aforementioned quenched investigations [8, 9] it was found that, e.g., the result 
for 6a — 6p at a ~ 0.1 fm does depend significantly on the computational details, namely the 
timeslice location xq and the exact lattice representation of derivatives. (See also the upper 
right panel in figure 2 of the present work.) Obviously, this is quite an extreme example of 
an improvement coefficient that is difficult to determine, and its ambiguity from choosing 
a different improvement condition may even grow to the same order as the coefficient itself. 
The reason for such a behaviour lies in the fact that the 0(a) effects are not dominating 
over the O(a^) ones in the considered correlation functions. Presumably, the 0(a) effects 
are just not very large. 

One of the benefits of the Schrodinger functional pointed out before is that improve- 
ment conditions formulated through SF correlation functions can also be studied in per- 
turbation theory. Choosing kinematical variables such as xq/L, T/L and 9 exactly as 
in the non-perturbative setup, the expansion of an improvement coefficient, here denoted 
generically by 6, is calculated as a series 

6(5o,a/i)~6W(a/L) + g2&('^(a/L) + 0(5o') . (2.32) 

For our example of 6a — 6p it was observed in ref. [8] that the perturbative coefficients 6^*^^ 
and 6*^^) show a similar dependence on the kinematics as the non-perturbative results. Our 
statement on the uniqueness of the improvement coefficients in perturbation theory at the 
beginning of this section, however, means that although generically the functions b^^\a/L) 
and b^^\a/L) will depend on the kinematical choices made in the improvement condition, 
the values 6'^*^)(0) and 6(^^(0) do not and are unique. 
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Hence, the following generic picture emerges. Once improvement conditions at constant 
physics are imposed, by which we mean that all length scales in the correlation functions 
that define the improvement condition are kept fixed in physical units, only the lattice 
spacing a changes when go is varied. This is just the situation, to which Symanzik's local 
effective theory of cutoff effects can be applied. As a consequence, improvement coefficients 
extracted from different constant physics conditions are then expected to rapidly approach 
an almost unique function of go as go ^ 0. While it is clear from this discussion that 
fixed L in physical units is to be preferred whenever possible, in certain situations one may 
want to implement improvement conditions at fixed L/a for practical reasons. But in such 
a setting, the 0(a) ambiguities in the improvement coefficients would not go to zero as 
90^0 (a^0).3 

In the present situation, where we have to work at finite quark masses, also the com- 
bination Lmn^i should be kept constant. Note, however, that these constant physics condi- 
tions do not have to be satisfied very precisely, because this would affect corrections to an 
0(a) term only. In addition to the studies [8, 9] in quenched QCD, such conditions have 
meanwhile also been imposed for the determination of ca [20] and Za [26] in the dynamical 
two-flavour theory. 

Let us summarize the discussion of this subsection. While an unfortunately chosen 
improvement condition, but implemented at constant physics, may even enlarge the cutoff 
effects for intermediate a, it guarantees that linear a-effects are absent. Putting the im- 
provement coefficient to zero or some perturbative approximation does not guarantee the 
latter, and the linear a-effects should at least be estimated in some way. 

3. Numerical computation and results 

The goal of this section is to non-perturbatively compute the improvement factors 6m (so) 
and [6a — bp]{go) as well as the renormalization constant Z(go) in the N{ = 2 dynamical 
case along the lines of the foregoing section. 

3.1 The constant physics condition 

In view of the freedom of possible legitimate choices for kinematical parameters within the 
SF formulation of lattice QCD [10], it was found to be advantageous to choose slightly 
different, particularly adapted settings in specifying a suitable constant physics condition 
for each of the parameter regions mentioned in the introduction. 

In the SF, the renormalization scale fi = 1/L and thus the physical volume are implic- 
itly defined by demanding the non-perturbatively renormalized coupling to take a certain 
value. Therefore, to impose our improvement/normalization conditions on a line of con- 
stant physics, it is required to keep all length scales fixed as go is varied, which we realize 
for the purposes of this work as 

^E.g., improvement conditions at fixed L/a were employed for tiie non-perturbative determination of 
Csw and Ca in quenched QCD [25]. In tfiis case, however, tlie small tree-level a-effect was subtracted from 
the non-perturbative ones to insure that the improvement coefficients go to their tree-level values exactly; 
furthermore, the conditions were chosen such that also \b^^\a/L) — fe'^'(0)| is very small. 
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• condition 1: 

(9 = 0.5, g'^{L) = 2.989 ^ L 0.25 fm for 6.2 < /3 < 6.6 (3.1) 

• condition 2: 

9 = 0, L^l.SL* s.th. g'^{L*) = 5.5 ^ L 1.2 fm for 5.2 < /3 < 5.7 (3.2) 

supplemented by choosing as fixed kinematical parameters 

Ck = C'j^ = (vanishing background field) , T/L = 3/2 , (3.3) 

where Ck and C'f, induce the Dirichlet boundary conditions of the gauge field at xq = 
and T, respectively. 

The first /3-range ('condition 1') refers to a region of weak couplings encountered in 
solving renormalization problems of lattice HQET and in studies of B-physics observables 
in the continuum limit of finite-volume two-flavour lattice QCD, see refs. [27, 28] and 
section 4. Our second /3-range ('condition 2') corresponds to lattice spacings of 0.09 fm > 
a > 0.04 fm [29, 30], which is the typical coupling region employed in A^f = 2 large- volume 
simulations. In this case, the choice 9 = (periodic boundary conditions) together with an 
intermediate physical volume of L ~ 1.2fm is guided by experiences made in the context 
of the non-perturbative redetermination of the axial current renormalization constant Za 
in ref. [26]. There it was observed that this is a kinematical region, where a^-effects 
(which can be significant at a ~ 0.1 fm) are better suppressed, implying reduced intrinsic 
a^-ambiguities in Za, and SF simulations can be done close to or in the chiral limit of 
vanishing sea quark mass. The chiral limit can be reached in the simulations, when the 
infrared cutoff on the spectrum of the Dirac operator (in the SF dominantly controlled by 
1/r) is sufficiently large, so that T = 3L/2 ~ 1.8 fm is a feasible choice for this parameter 
region and the underlying A^f = 2 gauge field configurations could be shared among the 
two projects. 

For the generation of our dynamical gauge field configurations with N{ = 2 massless 
sea quarks, we employed the HMC algorithm with two pseudo-fermion fields as proposed 
in [31, 32]. The particular implementation has been discussed and tested in refs. [33, 34] 
and following the latter, for the /3-range and volume of (3.2) we always chose a molecular 
dynamics trajectory length of r = 2, the only exception being /? = 5.2, where we set r = 1. 
For the /3-range and volume of (3.1), we had r = 1 throughout. To prepare a light sea 
quark doublet of mass zero in the simulations, the corresponding hopping parameter, ki, 
was tuned to its critical value, Kc- 

In order to complete the definition of the lines of constant physics, values for the bare 
PCAC masses of the valence quarks must be selected. In order to get an impression of the 
quark mass dependence in physical observables, we consider in analogy to the quenched 
study [9] two sets, 

set 1 : Lniii « , Lm22 ~ 0.5 , (3.4) 

set 2 : Lniii w , Lm22 ~ 2.5 , (3.5) 
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LI a 


Tja 


/3 






Lmu 


10 


15 


6.1906 




0.1360160 


-0.0054(7) 








K2 


0.1343182 


+0.5004(6) 








K2 


0.1276218 


+2.5004(7) 


12 


18 


6.3158 


Kl 


0.1357930 


-0.0062(10) 








K2 


0.1343782 


+0.5008(9) 








K2 


0.1287549 


+2.5007(10) 


16 


24 


6.5113 


Kl 


0.1354410 


-0.0101(13) 








K2 


0.1343872 


+0.4949(13) 








K2 


0.1301456 


+2.4955(13) 


20 


30 


6.6380 


Kl 


0.1351630 


+0.0176(9) 








K2 


0.1343562 


+0.5024(9) 








K2 


0.1309650 


+2.5039(10) 



Table 1: Results after tuning the hopping parameters in our smah- volume computations to different 
values of the dimensionless (and plateau-averaged) PCAC masses Lmx\ and Lm-ii as specified in 
eqs. (3.4) and (3.5). Set 1 and set 2 share the same sea quark content, given by k\ = Kc = The 
hopping parameters marked by a hat belong to set 2. 

where throughout this section all bare current quark masses are understood to be calcu- 
lated according to an averaging prescription analogous to eq. (2.31). The appearing local 
current quark masses raaix^) are evaluated using the standard lattice derivative as implied, 
following standard conventions, by its definition (2.20) based on the 0(a) improved PCAC 
relation (2.9). 

The choice on Lm22 in set 1 is motivated by experiences from the quenched investi- 
gation, ref. [8], where it was argued to be advantageous w.r.t. the size of 0(a) ambiguities 
encountered in perturbation theory. By contrast, the choice of Lm22 ~ 2.5 in set 2 is 
closer to the typical b-quark region itself. The conditions on Lm22 have been satisfied 
by adjusting the hopping parameters K2, responsible for the mass value of the heavy va- 
lence quark flavour. The hopping parameter needed to extract 6m, which fixes a third 
mass value according to eq. (2.30), is given by K3 = 2k,ik,2/{ki + K2) and was directly 
implemented in the measurement code to avoid rounding errors. Furthermore, we exploit 
the additional time-reflection symmetry of the forward and backward boundary-to-bulk 
fermionic SF correlation functions introduced in the previous section by summing them in 
the analysis accordingly. 

3.2 Results in the w^eak-coupling region of finite-volume QCD 

The parameter values for calculating the ratios Rx in L ~ 0.25 fm are given in table 1. 
Our statistics consists of 200 — 300 gauge field configurations stored at frequencies 5 — 10 
in units of molecular dynamics trajectories. The time extent T = 3L/2 is chosen such 
that the plateau average, eq. (2.31), is less affected by finite-size effects stemming from 
the boundaries compared to the other natural choice T = L. After having produced gauge 
field configurations with Nf = 2 massless dynamical quarks, we performed one or two steps 
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Figure 1: Time dependence of dimensionless PC AC masses obtained in small- volume QCD simula- 
tion after tuning the hopping parameters of the valence quarks as given in table 1. The dashed lines 
define the plateau region used for the averaging as described in the text. Corresponding averaged 
mass values are shown as error bands in each subplot. Lmn « 0, and Lm22 ~ 0.5 for set 1 resp. 
Lm22 ~ 2.5 for set 2 are fulfilled within the desired accuracy. 



of tuning to properly adjust K2 and satisfy the conditions'^ imposed on Lm22 (see also 
figure 1). The estimators Rx, X G {AP,m, Z}, are defined in terms of PCAC masses, 
which are secondary quantities themselves, built from SF correlation functions by means 
of eq. (2.20). For the evaluation of all our observables and their statistical errors we hence 
employed the F-method [35], which directly analyzes their autocorrelation functions. As a 
cross-check, a standard Jackknife procedure was applied as well. The results for Rx from 
improvement condition 1 for set 1 and set 2 are listed in table 2. 

The two plots on the r.h.s. of figure 2 show the timeslice dependence i?Ap(a^o) and 
Rm{xQ) for set 1 at L/a = 12 and /3 = 6.3158. The averaged estimators that represent 
our main results are shown as error bands. The filled symbols refer to improved time 
derivatives, (2.21), while the open symbols are obtained using the standard derivatives. 
Their difference can become quite large especially for Rap-, as can be inferred from the plot. 
This is expected to be a pure cutoff effect and underlines the necessity to use improved 
derivatives. A plateau in Rx^{xq) is visible in the data for both, standard and improved 
derivatives, but it is more pronounced for the latter. 



^Similar to the situation in [8, 9], this is to a sufficient precision equivalent to keeping the corresponding 
renormalized masses LZpjnjZ-p fixed, as for the considered couplings the entering overall renormalization 
constant barely varies. 
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set 


L/a 


P 




fern 


Z 




1 


10 


6.1906 


+0.0027(11) 


-0.6576(15) 


+1.10418(19) 


+0.6637(8) 




12 


6.3158 


+0.0011(37) 


-0.6666(26) 


+1.10522(29) 


+0.6653(13) 




16 


6.5113 


-0.006(11) 


-0.6637(36) 


+1.10395(29) 


+0.6614(22) 




20 


6.6380 


-0.0050(29) 


-0.6636(48) 


+1.10333(23) 


+0.6648(21) 


2 


10 


6.1906 


+0.07261(41) 


-0.56417(38) 


+1.09522(11) 


+0.63820(37) 




12 


6.3158 


+0.05177(61) 


-0.57800(66) 


+1.09747(13) 


+0.62971(60) 




16 


6.5113 


+0.02950(60) 


-0.5955(10) 


+1.09945(13) 


+0.62617(69) 




20 


6.6380 


+0.02101(68) 


-0.6077(11) 


+1.09997(11) 


+0.62914(71) 



Table 2: Summary of results on the improvement coefficients and tfie renormalization factor under 
consideration in the volume with L « 0.25 fm, for both sets of heavy quark masses (resp. constant 
physics conditions) using plateau-averaged masses and improved derivatives. 

For a careful error analysis we also looked at the dependence of the mean values and 
errors on the number of points entering the plateau average about xq = T/2. This is 
exemplarily shown in the left plot of figure 2. There the leftmost point is the one at 
xo = T/2 in the plots on the right. In the examples given here, the plateau average is 
approached fast. This also nicely reflects that a plateau average is less affected by local 
fluctuations, which could distort the mean value by more than its associated uncertainty. If 
a plateau is visible, the averaging prescription — being considered as part of the definition 
of the Rx as pointed out before — leads to a more reliable estimate for the mean value 
and the corresponding error. 

The dependence of our main estimates Rx is visualized in the left plot of figure 3. The 
non-perturbative data points show a smooth dependence on as expected, suggesting that 
the correct scaling of discretization errors proportional to a for the improvement coefficients 
and a? for the renormalization constant is met. Red data points correspond to our preferred 
choice, set 1, whereas the blue points belong to set 2. The qualitative dependence on Qq is 
comparable to the quenched case [9]. Note, however, that such a direct comparison at the 
quantitative level is of limited meaning only, since in the latter case the bare parameter 
region of the simulations referred to a smaller physical volume and the relation between 
the lattice spacing and the bare coupling is different from the two-flavour theory. 

In our fixed-volume computations, the lattice resolution gets smaller with decreasing 
Qq. Accordingly, the difference between both data sets decreases in this direction as it 
should be. Except for i^APi set 1, the non-perturbative points show a partly significant 
deviation from the one-loop perturbative predictions [8, 36] that are drawn as dashed lines 
in figure 3, too. For R^ this deviation is of 0(10%) and for Rz still about 1%. One 
also observes the results of set 2 to have a larger curvature, but a smaller error in the 
range of bare couplings considered here. This may be explained by the much larger valence 
quark mass, Lm22 \ set2 ~ 5Lm22 | seti, where the dynamics is governed by the heavy quark 
flavour and dominates the fluctuations induced by the dynamical massless quarks in the 
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Figure 2: A typical data sample, showing Rap and i?m data at L/a = 12 and f3 = 6.3158 for set 1. 
Left: Dependence on the plateau size used for the averaging, with the preferred choice indicated by 
dotted hnes. Right: Timeslice dependence of Rx{xo), if no plateau averaging of the local PCAC 
masses is applied. Open symbols denote results using standard derivatives, while filled symbols 
refer to those with improved lattice derivatives, eq. (2.21). For comparison, the shown (symmetric) 
error bands correspond to the results obtained from plateau-averaged masses. 



gauge background. 

The difference of Rx between two improvement conditions is a measure for its 0(a") 
ambiguity. In the extreme case of bA — bp, this difference in the results between sets 1 
and 2 is of the same order as the improvement coefficient itself so that in such a situation 
it appears rather tempting to just put the coefficient to zero. In the light of the qualitative 
discussion in subsection 2.3, however, it should be obvious that it is then not guaranteed 
that linear a-effects are absent after improvement. Thus, an extrapolation to the continuum 
using an 0{a?) model function for the cutoff effects might then give significantly wrong 
results. Based on this warning, let us now study such ambiguities quantitatively. 
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Figure 3: gQ-dependence of our final results for 6a ~ &p and from top to bottom. The dashed 
lines refer to the predictions from one-loop perturbation theory. Left: Results in the weak-coupling 
region of small-volume QCD (constant physics condition 1) for set 1 (red points) and set 2 (blue 
points). Right: Same points together with the corresponding results of subsection 3.3 in the region 
of larger couplings commonly employed in large-volume simulations (constant physics condition 2, 
green points). 

0(a") ambiguities 

As a consequence of the very idea of the constant physics condition, any other estimate 
Rx (i.e. stemming from a different choice of improvement or renormalization condition) in 
general will yield a different functional dependence upon (and therefore a/L), but its 
differences are again smooth functions that must vanish in the continuum limit with a rate 
proportional to a/L for improvement coefficients or (a/L)'^ for renormalization constants. 
Regardless of their actual size, these intrinsic 0(a") ambiguities (n = 1,2) imply that 
rather than a numerical value at some given f3, the important information lies in the 
correct gQ-dependence of the estimators Rx, X £ {AP,m, Z}, obtained by working at 
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Figure 4: The differences ARx of the estimators Rx among the two sets of results, as defined 
in eq. (3.6). Within errors, the differences approach zero with a rate proportional to a/L (for 
improvement coefhcients) and {a/LY (for the rcnormalization factor) as expected. The earlier 
quenched results [9] are added for comparison. 



constant physics while varying j3. To demonstrate this, we investigated a few alternative 
improvement conditions, which are either already provided by the two quark mass choices 
themselves, eqs. (3.4) and (3.5), or by defining the estimators Rx with standard instead 
of improved derivatives. 

Regarding the latter, we mentioned before that cutoff effects are visible between the 
use of standard and improved time derivatives, which are used to build the PCAC masses, 
and the reader may have wondered already about the large differences between these two 
choices of lattice derivatives in figure 2. But one should keep in mind that these results 
refer to an intermediate lattice spacing with L/a = 12. When performing the full analysis 
for each data set with the standard derivatives, however, the differences to the results from 
employing improved derivatives were found to vanish in the limit a — )■ as expected. 

In addition, the two choices (3.4) and (3.5) made for the non-vanishing valence quark 
mass themselves correspond to different improvement conditions to extract the estimators 
Rx and so may yield an impression of the intrinsic cutoff ambiguities. To this end we 
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Figure 5: Unconstrained continuum extrapolations of ii/ps(2;, .go) for two fixed values of the 
renormalized quark, with different and Z entering the latter as described in the text. The data 
obtained via set 1 (red squares) and via set 2 (blue triangles) is fitted using three points (solid lines) 
and four points (dashed lines), respectively. 



From eqs. (2.26) and (2.27) one expects a linear a-dependence and according to (2.28), 
the cutoff dependence of Rz should be of O(a^). For the values listed in table 2 we show 
the results for ARx as red points (diamonds) in figure 4. While the difference ARz in 
the right panel clearly exhibits a linear approach towards zero, the 0(a) ambiguities of 
6a — bp and reflected in the left panel are quite small as well, and their magnitudes 
rapidly decrease as a/L — )• 0. The additional black points (circles) obtained in the quenched 
case [9] lie in the same ballpark, which can be attributed to the fact that the impact of the 
dynamical light quarks largely cancels out in the differences (3.6); on top of that, the cutoff 
ambiguities are dominated by the heavy quark mass content, which is about the same in 
both studies. In summary, the scaling behaviour observed in the various differences of 
estimators (3.6) confirms the theoretical expectation of the universality of the continuum 
limit and supports the reliability of the applied non-perturbative method also in the case 
of A^f = 2 dynamical quarks. 

Although based on theoretical grounds discussed in subsection 2.3 we expect that phys- 
ical observables involving sets of improvement coefficients from different constant physics 
conditions will have the same continuum limit, let us illustrate this in a concrete example. 
We consider the dimensionless combination Li/ps, deriving from an effective energy that 
is extracted from the (improved) SF correlator of the time component of the axial vector 
current in a finite volume Lf with Li ~ 0.5 fm, 



- ..... ^^^^^^2 ^ V , 

We confront its lattice spacing dependence using and Z of set 1 with that using and 
Z of set 2 to fix the dimensionless renormalization group invariant (RGI) heavy quark mass 



define 





(3.6) 




z = LiM . 



(3.7) 



to z = 12, 16 for lattice resolutions of Li/a = 20 — 40. For further technical details, see 
section 4. Both data sets and its unconstrained continuum extrapolations linear in (a/L)^ 
are displayed in figure 5 and extrapolate to the same continuum limit. Associated statistical 
errors are smaller or of the order of the symbol size. The assumed quadratic scaling 
behaviour is clearly violated for the coarsest lattice at z = 16 in set 1 though, because 
quark mass dependent cutoff effects are substantial in that case. Hence, we discarded 
the L/a = 20 points in our fits of data set 1. The fact that cutoff effects in Lilps for 
both cases are larger for set 1 is not surprising, since the improvement condition of set 2 
with Lom22 ~ 2.5 and thus Lim22 ~ 5 is much closer to the line in parameter space 
with z = 12, 16 along which Li/ps was computed. This general behaviour present in our 
new data was also observed in the quenched case [9]. The slope of the blue solid line at 
z = 16 is approximately 6.5 times larger than the corresponding one at z = 12. However, 
from the excellent agreement of the continuum limits in both cases one infers once more 
that our results correctly model the ^Q-dependence of and Z. No matter which set 
of improvement conditions is used, both entail convergence to the continuum limit with 
leading corrections of O(a^) in the 0(a) improved theory. Consequently we can conclude 
that the ambiguity introduced by choosing either set 1 or set 2 vanishes in the continuum 
limit also for other values of z and/or observables and universality holds. 

3.3 Results in the coupling region of large-volume simulations 

Besides in the parameter range discussed in the previous subsection, which is relevant, for 
instance, in the context of the matching calculations of HQET and QCD and to which we 
will come back again in the next subsection, the estimators Rx, X € {AP,m, Z}, are also 
of strong interest in the parameter region, where lattice computations of hadron masses 
and matrix elements in physically large volumes usually take place. This is particularly the 
case for mesons involving the charm quark, because large-volume studies in the quenched 
approximation [17, 18] have already given clear evidence for possibly substantial quark 
mass dependent cutoff effects down to a ~ 0.05 fm so that a non-perturbative removal of 
the associated 0(a) corrections is highly desirable. Thus, we here provide 6a — 6p, and 
Z as functions in covering the range 5.2 < /3 < 5.7, which then will allows to interpolate 
to the typical /3- values relevant for large- volume computations in two-flavour lattice QCD. 

In table 3 we summarize the corresponding simulation parameters, where the relation 
between L* /a and from ref. [30] was taken to translate the constant physics condition 2, 
eq. (3.2), to pairs of {L/a, (3). Owing to algorithmic instabilities caused by the appearance 
of very small, unphysical eigenvalues of the SF Wilson-Dirac operator [37, 38] , the lattice 
with the coarsest lattice spacing (i.e. L/a = 12 at /3 = 5.2) could only be simulated down to 
bare PCAC sea quark masses of about am\ = am^^°^^ ~ 0.01. Therefore, our computations 
for this parameter set were based on dynamical configurations at three sea quark mass 
values > 0.01, made available to us through the determination of Z\ reported in [26], in 
order to extrapolate the results for /3 = 5.2 to the chiral limit. 

To eventually extract the estimators Rx through eq. (3.2), we here considered the 
setting (3.4) for the non-degenerate valence quarks so that the heavier PCAC quark mass 
Lm22 had to be tuned to ~ 0.5, while the light one was again identified with the sea quark 
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L/a T/a (3 ki Lmn K2 Lm22 total 



12 



18 5.20 



0.135500 
0.135700 
0.135800 



0.2495(30) 
0.1584(32) 
0.1051(30) 




0.134677 
0.134601 
0.134512 



0.5006(28) 4000 

0.4956(28) 4000 

0.5015(23) 4000 
04986(36) 



16 24 5.40 0.13664 5 0.0094(26 ) 0.13554 ? 0.4996(20) 750 
24 36 5.70 0.136704 0.0141(10) 0.136000 0.5082(9) 656 

Table 3: Simulation parameters to compute and 6a — &p and Z for a range of bare couplings 
typically covered in large-volume simulations of A'f = 2 lattice QCD. The rightmost column gives 
the number of decorrelated gauge field configurations used. The kinematical SF parameters are 
(T/L, 61) = (3/2,0), cf. eqs. (3.2) and (3.3). 



r/3 avr'd xq = T/2 r/3 avr'd 



T/2 r/3 avr'd xq = T/2 



Ki 

Lrriii 
TTO22 

Z 



0.135800 
0.1051(30) 0.1072(47) 



0.5015(23) 

-0.098(50) 
-0.453(28) 
0.9262(23) 



0.5015(37) 

-0.078(95) 
-0.413(43) 
0.9194(38) 



0.135700 
0.1584(32) 0.1668(47) 



0.4956(28) 

-0.040(42) 
-0.513(25) 
0.9305(21) 



0.5016(39) 

-0.048(74) 
-0.511(43) 
0.9225(35) 



0.135500 
0.2495(30) 0.2547(43) 



0.5006(28) 

-0.086(25) 
-0.485(20) 
0.9191(16) 



0.5041(37) 

-0.108(45) 
-0.424(36) 
0.9058(31) 



Table 4: Results at /3 = 5.2 before extrapolating to the massless limit, stemming from averages 
over the central T/3 timcslices and from using .tq = T/2 only. 



mass, Lmii = Lm\, by setting ki = k\ (= Kc in case of /3 = 5.4, 5.7 only). Lmu and Lm22 
are also listed in table 3. 

For the (3 = 5.2 data we do not observe a pronounced plateau in Rx(xo), and we thus 
gather in table 4 results stemming from our standard definition with an average taken over 
timeslices of the central third of the time extent, together with the corresponding results 
obtained at xq = T/2. Note that the improvement coefficients are well compatible within 
errors, whereas for the renormalization constant this is only the case for the smallest light 
quark mass. We hence decided to not alter our procedure of estimating these quantities 
using timeslice averages. As for the chiral extrapolation Lmu = Lmi — ?• of the /3 = 5.2 
results, we investigated fits of Rx to a constant as well as linear and quadratic fits using 
two or all three (light valence = sea) quark mass points, which gave overall consistent 
results, since the dependence on the light quark mass is rather mild. Two- and three-point 
weighted averages turned out to be most conservative and as estimates in the massless 
limit we just quote the weighted averages of the numbers at the two lightest quark masses. 
To account for a systematic error due to O(ami) contaminations, we follow ref. [26] and 
inflate the error of the mean by its difference to the heaviest sea quark mass point linearly. 
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Figure 6: Extrapolations of the results at /3 = 5.2 in table 4 to the chiral limit, Lmu = Lm\ = 0. 
The weighted two-point averages shown in the plot include our estimate of the 0{am\) systematic 
uncertainty as explained in the text and quoted in table 5. 



L/a (3 



6a - bp 



Z 



12 5.20 0.135500 
0.135700 
0.135800 



-0.086(25) 
-0.040(42) 
-0.098(50) 
-0.064{32)i22) 



-0.485(20) 
-0.513(25) 
-0.453(28) 
-0486{19){1] 



0.9191(16) 
0.9305(21) 
0.9262(23) 
0.9285{16)i94) 



16 5.40 0.136645 -0.159(37) -0.562(36) 1.0319(22) 
24 5.70 0.136704 - 0.014(21) -0.770(21) 1.0976(8) 

Table 5: Results of 6a ^ and Z in the bare parameter region 5.2 < /3 < 5.7. The quoted 

error on the results for /3 = 5.2 is split into a statistical and a systematic one, the latter accounting 
for possible 0(a77ii) contaminations from the extrapolation to the massless limit. 



Our final results for all values of (3 are summarized in table 5. Since simulations with 
massless dynamical quarks in this parameter region are numerically much more expensive, 
the P = 5.4 and 5.7 data samples have a much smaller statistics compared to the ones 
reported in subsection 3.2. Consequently, the estimates in this parameter region have an 
uncertainty, which in general is one order of magnitude larger than in our small volume 
simulations at weaker couplings. Even though one might generically expect a universal 
curve connecting the results in distinct parameter regions, one encounters significantly 
different cutoff effects owing to different improvement conditions applied, as reflected on 
the r.h.s of figure 3. 

For later use of our results it seems appropriate to represent them as smooth func- 
tions in the interval 5.2 < /3 < 5.7. As in the quenched computation [8], we can readily 
parametrize Z, — bp and 6m by means of Pade fits in Qq that are constrained to pass 
asymptotically into the perturbative one- loop predictions in the limit Qq ^ 0. In case of 
the renormalization factor Z, the fit 

Z (/„) = (1 + 0,0905149?) X I -O -^^^^oi- 0.214^4 p.g, 
^ 1- 0.6186 c/^ ^ ^ 

is found to describe our numbers very well, and the error to be attributed to this formula 
varies from about 1.2% at /3 = 5.2 to less than 0.1% at /3 = 5.7. Due to the scatter in the 
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results for 6a — bp, we quote as Pade fit parametrization 



[bA-bp] (gl) = -0.00093 x ^ + '^^■^^9^ (3.9) 
I A fi yyoj yo 1 _ 0.6235 5^ ^ ^ 

which represents the data within two sigma of their overall uncertainty. Note that acci- 
dentally the parameter in the numerator is close to the quenched value, 23.306, and yet an 
unconstrained linear fit would represent the data points nearly equally well. Finally, for 
the improvement coefficient bm we quote the interpolating formula 

I — Q 3737 

bra (50) = (-0.5 - 0.09623 5o') x ^_^^^^^gA ' (3-10) 

still respecting the leading perturbative asymptotics. When using this formula at i^Q-values 
different from (3 = 5.2,5.4,5.7, one should attribute to an uncertainty of the order of 
the error quoted in table 5 for the result at the nearest /3. 



4. Lines of constant physics in finite-volume QCD 

As an immediate application of the results obtained in subsection 3.2, we address the 
computation of the dependence of heavy-light meson observables on the heavy quark mass 
in the continuum limit of finite-volume QCD. This is one of the key elements in the basic 
strategy [14] to non-perturbatively determine HQET parameters by a matching to QCD. 
Here, we only restrict the discussion to the fixing of the RGI heavy quark mass in a given 
finite volume to a series of values covering the b-quark mass down to the onset of the charm 
quark region, because it is this part, where and Z directly enter. A preliminary account 
on the subsequent computation of the heavy quark mass dependence of finite- volume QCD 
observables was given in the report [28]. 

Since we are eventually interested in quark mass dependent quantities evaluated in 
the continuum limit of finite- volume QCD, we want to calculate them for a series of bare 
parameters (L/a, /3, k\) such that the renormalized parameters in the light quark sector are 
fixed and in this way physics is kept constant along the approach to the continuum limit. 
As before, k\ is the hopping parameter of the mass degenerate dynamical light quarks, 
while their PCAC mass is denoted by mi. In order to also fix the renormalized masses 
of the non-degenerate valence quarks in units of a physical scale later on, one first has to 
specify a suitable value of the renormalized SF coupling, g'^{L), which then, by its very 
definition [10], automatically fixes the spatial extent of the SF cylinder as the only scale of 
our system. Thus, in addition to 9 = 0.5 and T = L, we choose 

g\Lo) = 2.989 , ^0 = T ' ""^^^^^ " ° ^^'^^ 

and thereby impose a constant physics condition on the renormalized coupling and the 
mass of the sea quark doublet in a finite volume of extent Lq. Implicitly, this condition 
defines the volume characterized by Li = 2Lq, in which the non-perturbative matching 
between HQET and QCD will be performed [27]. Although an exact knowledge of Li in 
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L/a 


P 


^p(ffo,i/a) 


set 




byn 


Z 


10 


6.1569 


0.6065(9) 


1 


-0.0000(12) 


-0.6633(12) 


1.10443(17) 








2 


+0.07852(53) 


-0.56196(38) 


1.09488(13) 


12 


6.2483 


0.5995(8) 


1 


-0.0016(8) 


-0.6661(9) 


1.10475(12) 








2 


+0.06284(33) 


-0.57145(29) 


1.09632(8) 


16 


6.4574 


0.5941(10) 


1 


-0.0050(17) 


-0.6674(23) 


1.10455(17) 








2 


+0.03567(51) 


-0.59147(63) 


1.09888(10) 


20 


6.6380 


0.5949(12) 


1 


-0.0045(28) 


-0.6692(27) 


1.10379(17) 








2 


+0.02150(58) 


-0.60763(78) 


1.10021(10) 



Table 6: Summary of the improvement coefficients and renormalization constants slightly extra- 
/interpolated to ensure a fixed SF coupling in a volume with L = Lg ~ 0.25 fm. The results 
at L/a = 20 had not to be extrapolated. Zp was determined independently following the non- 
perturbative renormalization procedure of [40]. 

physical units is not yet needed at this stage, one can infer from the known running of the 
SF coupling for iVf = 2, c.f. refs. [39, 30], that Li = 2Lo ^ 0.5 fm. Obviously, eq. (4.1) is 
identical to the corresponding settings of the constant physics 'condition 1', eq. (3.1), used 
for the determination of — bp, 6m and Z themselves in subsection 3.2 and, in retrospect, 
explains this so far supposedly arbitrary choice. 

The peculiar value ^^(Lq) = 2.989 was fixed by a simulation at Lo/a = 20, T = Lq, in 
conjunction with tentative interpolations in f3 = Q/g^ for given L^/a < 16 to this target, 
based on the known dependence of the SF coupling and the current quark mass on the bare 
parameters (/?, k) available from the data of ref. [39].^ The condition of zero light sea quark 
mass in eq. (4.1) is met by setting ki to the critical hopping parameter, Kc, estimated again 
on basis of published data [39], whereby a slight mismatch of \LQm\(LQ)\ < 0.05 of this 
condition is tolerable in practice. The triples (Lg/a, /3, ki) resulting from this procedure are 
those of table 1, used in our study of improvement and renormalization factors described 
above. In order to avoid a non-negligible systematic error from small violations of the 
condition (4.1) on results of any subsequent computations in Lf though, the /3- values were 
checked and fine-tuned further by additional simulations, which at the end led to the values 
in the second column of table 6. The corresponding small shifts in /? affect our previous 
estimates of bp^ — bp, ftm and Z only at a negligible level, because any deviation from the 
line of constant physics just entails a small change of the O(a^) effects that we nevertheless 
account for by interpolations in /3. The resulting, slightly moved values of b^ — bp, bm and 
Z can be found in table 6, too. 

For the intended study of finite-volume QCD observables [28] we consider SF corre- 
lation functions composed of a non-degenerate heavy- light valence quark doublet, where 
the light valence quark mass is again chosen to be equal to the vanishing sea quark mass 

^Using the known /^-function and our experience from the quenclied calculation [41], we can estimate 
that an uncertainty of about 0.04 in the couphng would translate via Li into an uncertainty in the b-quark 
mass of at most 0.5%. 
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adopted in the generation of the dynamical gauge configurations, i.e. it is set to ki = Re- 
computing the heavy quark mass dependence of observables then amounts to evaluate 
the correlators on these configurations in a volume Lf x T, T = Li, for a sequence of 
renormalized heavy valence quark mass values in a region extending from charm to beauty. 
For this purpose, the latter have to be fixed as precisely as possible for the same given 
pairs (Lo/a, /3) that already fix the volume Li via enforcing constant physics aX, Lq = Li/2 
through eq. (4.1). Therefore, starting from (2.13), we exploit the 0(a) improved relation 
between the subtracted bare heavy quark mass mq^h and the corresponding RGI mass, viz. 



M = h{Lo) Zraigo, Lo/a) (1 + 6m(9o) amq^h) "iq,h + O(a^) , (4.2) 



where 



7 r T / \ ^(90) ZA{go) 1/1 1 \ o\ 
Zrn{go,Lo/a) = — — ^ , amq^h = , (4.3) 

and the axial current renormalization constant Za in dependence of 50 is non-perturbatively 
know for = 2 from ref. [26]. The scale dependent renormalization constant Z-p was non- 
perturbatively determined in the course of our simulations in L = Lq, following exactly the 
prescription in ref. [40]. Its values relevant here (modulo a slight extrapolation in /? as in 
the case of 6a ^ 6p, bm and Z) are listed in the third column of table 6. In eq. (4.2), there 
also appears the factor 

M 12 
HLq) = —— = 1.521(14) , = (4.4) 

m(/io) Lq Li 

which represents the universal, regularization independent ratio of the RGI heavy quark 
mass, M, to the running quark mass m in the SF scheme at the renormalization scale fiQ. 
h{Lo) was evaluated by a reanalysis of available data on the non-perturbative quark mass 
renormalization in two- flavour QCD, also published in ref. [40]. 

Having all ingredients of eqs. (4.2) - (4.4) at hand, particularly the improvement 
coefficient 6m and the renormalization constant Z obtained in this work, these equations 
can now straightforwardly be solved for any specific value of the dimensionless RGI heavy 
quark mass in Li, 

z = LiM , (4.5) 

to yield the desired hopping parameter associated with the mass of the heavy valence quark 
flavour, 

for which the heavy-light SF correlation functions and meson observables composed of 
them are to be evaluated. For 6m (90) and Z{gQ) we employ the numbers belonging to the 
'set 1' (cf. eq. (3.4)) of our results, because firstly, the ^Q'dependence is weaker for 6m and 
secondly, with the setting LQni22 ^ 0.5 in their determination (see eq. (3.4)) we rather meet 
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K4 


9 


0. 
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0. 
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0, 
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0. 
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K5 


11 


0. 
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0. 


1278628 


0, 


1298401 


0. 


1307738 


Kg 


13 


0. 


1234098 


0. 


1261106 


0, 


1286909 


0. 


1299065 


K7 


15 


0. 


1204339 


0. 


1241815 


0, 


1274876 


0. 


1290126 


Kg 


18 






0. 


1206988 


0, 


1255509 


0. 


1276125 


Kg 


21 






0. 


1140810 


0, 


1233865 


0. 


1261232 



Table 7: Hopping parameters Ki = from eq. (4.6), corresponding to fixed dimensionless RGI 
heavy quark masses z = LAI in a physical volume with L = Li ^ 0.5 fm for different lattice 
resolutions. 

the condition that in ref. [8] was found advantageous from the point of view of perturbation 
theory^. Owing to the sign and the order of magnitude of the non-perturbative values for 
ftrn in the /3-range relevant here, eq. (4.6) has no real solutions for arbitrarily high z-values. 
This demands to balance the upper bound on z against the lattice resolutions one wants 
to simulate. In our case of Li/a = 2Lo/a = 20, 24, 32, 40 we thus opted for nine values, 

z G {4, 6, 7, 9, 11, 13, 15, 18, 21} , (4.7) 

which fairly span a range from around the charm to beyond the b-quark mass region. 
Equation (4.6) then yields the corresponding sets of heavy-flavour hopping parameters for 
each Li/a collected in table 7, where solutions of (4.6) with z = 18, 21 only exist for inverse 
lattice spacings of Li/a = 24, 32, 40. 

Verification of 0(a) improvement 

For fixing z to the values in (4.7), one specific definition of the quark mass was applied, 
namely the relation between the bare subtracted heavy quark mass and the corresponding 
RGI mass (4.2), where the previously determined improvement coefficient bmido) and the 
renormalization constant Z{go) enter. Nevertheless it appears tempting to have a look at 
equivalent definitions of the RGI heavy quark mass based on the PCAC relation, in order to 
provide another check of universality of the continuum limit in the 0(a) improved theory. 

®We recall in this context that it was explicitly demonstrated in subsection 3.2 that a typical observable 
such as the effective heavy-light meson energy approaches the same continuum limit, when computed at a 
fixed value of z in which either bm and Z from 'set 1' or from 'set 2' enter. 
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0.002 0.001 0.001 0.002 

Figure 7 : Continuum extrapolations of the dimensionless RGI PCAC heavy quark masses, defined 
through either the heavy-light (left) or the heavy- heavy (right) axial current, in a volume with 
L = Li for 9 = 0.5. The statistical error is smaller than the symbol size. Dashed lines extend the 
fits to the points, which were not included in the extrapolation, and the bold ticks on the z-axis 
indicate the subset of z-values fixed via eq. (4.2) considered here. 



As a spin-off, this allows to also apply our non-perturbative knowledge of [b^ — bp] (go) from 
subsection 3.2 in an actual calculation, even though the influence of the results on 6a — bp 
belonging to the 'set 1' in this parameter region is expected to be numerically very small. ^ 
For each (Li/a, /3, ni), we take the heavy valence quark hopping parameters corre- 
sponding to the subset z S {4, 7, 11, 15, 18, 21} and compute the bare PCAC masses niij 
of eq. (2.20) from the heavy-light (i = h,j = \) and heavy-heavy (i=j = h) axial current 
and pseudo-scalar density. With eqs. (2.9) - (2.11) in mind, the renormalized PCAC quark 
masses then read 

5 ("T-R.i + m-nj) = -7^1 T—r\ {l + [6a - 6p](9o) 5 iam^^i + am^,j) }mij ■ (4.8) 

Since we are interested in recovering the above values of z = LiM by a continuum extrapo- 
lation of suitable PCAC quark masses, we need to express eq. (4.8) by RG invariants accord- 
ingly. Introducing the total renormalization factor Zyiigo) = h{LQ)ZA{gQ)/Zp{gQ, Lq/u), 
with h{Lo) from eq. (4.4), we end up with 

\{zi + Zj) = Zm{1 + [6a - bp] \ (arriq^i + am^j) }Limij . (4.9) 

From the PCAC quark mass definition via the heavy-heavy current one thus directly has 

z^"^^^ = Zm{1 + [6a - 6p] amq,h}Limhh , (4.10) 

while alternatively, z can also by obtained through the quark mass definition via the heavy- 
light current as 

z^^^') = 2Zm{1 + [6a - 6p] I (amq,h + am^^i) jLimu - z^^) ; (4.11) 

^Instead using the somewhat larger values of 'set 2' would not alter the outcome of this paragraph 
because A7?ap vanishes in the continuum Hmit. 
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here, z^^^^ as in eq. (4.10), but with h — t- 1. Since we have ki = Kc for our hght sea and 
valence quarks, this reduces to 

z^^^^ = 2Zm{1 + [6a - bp] iamq,h}^imhi . (4.12) 

Based on universaUty arguments, we expect the continuum extrapolations of eqs. (4.10) 
and (4.12) to approach the corresponding values of z fixed via the subtracted heavy quark 
mass, eq. (4.2). This is illustrated in figure 7. We show the results obtained from the 
heavy-light (z^^'^) and the heavy- heavy (z^^^^) current and compare the estimators for 
the PCAC masses using standard and improved lattice derivatives as introduced earlier, 
c.f. eq. (2.21). In qualitative agreement with our findings before, we in general observe 
that the cutoff effects increase with increasing values of z. The other expectation that 
cutoff effects are reduced when using improved lattice derivatives is clearly reflected as 
well. Whereas the continuum extrapolations for both lattice derivatives in the heavy-light 
case nicely match each other and coincide with the z-targets even at large values of z, one 
observes a slight mismatch for z > 15 in the heavy- heavy current case. This pattern can 
be explained by increasingly large quark mass dependent o-effects inherent to the standard 
lattice derivative. Since we used improved lattice derivatives in the course of estimating 
^A — bp, 6m and Z underlying this analysis, it does not come as a surprise that the smallest 
cutoff effects are seen for zC^^) with improved derivatives. This is further supported by 
the fact that the extrapolating fit function meets the data points even beyond the range 
included in the fit. 

Error budget 

To summarize, {Li/a, /?, ki, Kj), with Kj = and i = 1, . . . , 9, finally constitute the param- 
eters of the numerical valence quark propagator calculations, from which the heavy-light 
correlation functions and the heavy quark mass dependence of the observables constructed 
from them are to be extracted. If for some reason (say, because of yet having used tentative 
results for the quantities entering (4.6)) a subsequent calculation was based on a marginally 
different set of k- values, the associated small mismatch in having fixed z for each Li/a can 
easily be corrected by a small interpolation. 

At this point it still remains to quantify the uncertainty of z, which comes into play 
through the various factors and coefficients it depends on and which itself has to be prop- 
agated into any secondary quantity that is regarded as a function of z. With the estimate 
AZx/Za ~ 0.36% as quoted in [26] and AZp/Zp, AZ/Z as well as A6m/6m taken from 
table 6, ones obtains according to the standard rules of Gaussian error propagation an 
accumulated relative error on z in the range 0.38% < (Az/z) < 0.41%. This holds for all 
values of z in (4.7) and inverse lattice resolutions Li/a, but without the contribution from 
the continuum factor /i(Lo). The latter represents with Ah/h = 0.92% the dominating 
source of uncertainty in the total error budget of Az/z = 1.01%, when the fixing of z 
proceeds as described here. 
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5. Conclusions 



In this paper we have non-perturbatively determined the renormahzation constant and 
the improvement coefficients relating the renormahzed current and subtracted quark mass 
in 0(o) improved two- flavour lattice QCD. To solve this rather technical and difficult 
numerical problem, improvement and normalization conditions on lines of constant physics 
were imposed, which ensures a smooth dependence of these quantities on the bare coupling 
go and the correct scaling of discretization errors proportional to o (for the improvement 
coefficients) and (for the renormahzation constant), respectively. In particular, our 
results open the way to a precision study of B-physics on the lattice, whose feasibility 
has already been demonstrated within the very useful test-laboratory of the quenched 
approximation over the last years [15, 16, 41, 42]. 

An important element within these approaches is the separation of the scales char- 
acterizing heavy-light systems and the introduction of small physical volumes in order to 
perform QCD calculations with relativistic b-quarks. These simulations can only be per- 
formed by knowing the renormahzation constants and improvement coefficients computed 
in the present paper for the dynamical two-flavour theory. As explained in section 4, with 
the help of the latter the renormahzation group invariant heavy quark mass in finite- volume 
simulations of the 0(a) improved theory can accurately be fixed to values from the b- down 
to the c-quark region. Moreover, our results are also required when computing charmed 
observables by making simulations on the same volumes used to study light flavour physics. 

In the A'^f = 2 dynamical case, the programme of calculating the HQET parameters 
at 0(l/mb) non-perturbatively by a matching to QCD, which already employs part of 
the results of this work, is well advanced and progressing fast and will allow to determine 
several phenomenologically interesting quantities such as, for instance, the b-quark mass 
and B-meson decay constant. As a further interesting direction for applications of the 
results of our work let us mention the non-perturbative tests of the HQET expansion in 
the spirit of the quenched investigation [43] , which may greatly improve our confidence in 
the use of the effective theory. This will be presented in a forthcoming publication. 

Acknowledgments 

We are indebted to G. M. de Divitiis, M. Delia Morte, H. Meyer, R. Petronzio, H. Simma, 
R. Sommer and S. Takeda for fruitful discussions and their contributions at an early stage of 
this work. We also would like to thank R. Sommer for a critical reading of the manuscript. 
This work is part of the ALPHA Collaboration research programme. We thank NIC/DESY 
and INFN for allocating computer time on the APE computers to this project as well as 
the staff of the computer centers at Zeuthen and Rome for their support. We acknowledge 
partial support by the European Community through EU Contract No. MRTN-CT-2006- 
035482, "FLAVIAnet". P. F. and J. H. also acknowledge financial support by the Deutsche 
Forschungsgemeinschaft (DEC) under grant HE 4517/2-1. Some of the results presented 
here have already been published in [44]. 



-27- 



References 

[1] D. B. Kaplan, Chiral Symmetry and Lattice Fermions, 0912.2560. 

[2] R. Frezzotti and G. C. Rossi, C'hirally improving Wilson fermions. I: 0(a) improvement, J. 
High Energy Phys. 08 (2004) 007 [hep-lat/0306014]. 

[3] K. Symanzik, Continuum limit and improved action in lattice theories: 1. Principles and 4>'^ 
theory, Nucl. Phys. B226 (1983) 187. 

[4] K. Symanzik, Continuum limit and improved action in lattice theories: 2. 0(N) non-linear 
Sigma model in perturbation theory, Nucl. Phys. B226 (1983) 205. 

[5] B. Sheikholeslami and R. Wohlert, Improved continuum limit lattice action for QCD with 
Wilson fermions, Nucl. Phys. B259 (1985) 572. 

[6] M. Liischer and P. Weisz, On-shell improved lattice gauge theories, Commun. Math. Phys. 97 
(1985) 59. Erratum: ibid. 98 (1985) 433. 

[7] G. M. de Divitiis and R. Petronzio, Non-perturbative renormalization constants on the lattice 
from flavor non-singlet ward identities, Phys. Lett. B419 (1998) 311 [liep-lat/9710071]. 

[8] ALPHA Collaboration, M. Guagnelli, R. Petronzio, J. Rolf, S. Sint, R. Sommer and 

U. Wolff, Non-perturbative results for the coefficients 6ni and 5a — bp in 0(a) improved lattice 
QCD, Nucl. Phys. B595 (2001) 44 [liep-lat/0009021]. 

[9] ALPHA Collaboration, J. Heitgcr and J. Wennckers, Effective heavy-light meson energies in 
small-volume quenched QCD, J. High Energy Phys. 02 (2004) 064 [liep-lat/0312016]. 

[10] M. Liischer, R. Narayanan, P. Weisz and U. Wolff. The Schrddinger functional: A 
renormalizable probe for non-abelian gauge theories, Nucl. Phys. B384 (1992) 168 
[hep-lat/9207009]. 

[11] S. Sint, On the Schrddinger functional in QCD, Nucl. Phys. B421 (1994) 135 
[hep-lat/9312079]. 

[12] S. Sint, One-loop renormalization of the QCD Schrddinger functional, Nucl. Phys. B451 
(1995) 416 [hep-lat/9504005]. 

[13] K. Jansen et. al., Non-perturbative renormalization of lattice QCD at all scales, Phys. Lett. 
B372 (1996) 275 [hep-lat/9512009]. 

[14] ALPHA Collaboration, J. Heitger and R. Sommer, Non-perturbative heavy quark effective 
theory, J. High Energy Phys. 02 (2004) 022 [hep-lat/0310035]. 

[15] M. Guagnelli, F. Palombi, R. Petronzio and N. Tantalo, Fb and two scales problems in lattice 
QCD, Phys. Lett. B546 (2002) 237 [hep-lat/0206023]. 

[16] D. Guazzini, R. Sommer and N. Tantalo, Precision for B-meson matrix elements, J. High 
Energy Phys. 01 (2008) 076 [0710.2229]. 

[17] ALPHA Collaboration, J. Rolf and S. Sint, A precise determination of the charm quark's 
mass in quenched QCD, J. High Energy Phys. 12 (2002) 007 [hep-ph/0209255]. 

[18] ALPHA Collaboration, J. Heitger and A. Jiittner, Lattice cutoff effects for /d^ with 

improved Wilson fermions - a final lesson from the quenched case, J. High Energy Phys. 05 
(2009) 101 [0812.2200]. 



-28- 



[19] ALPHA Collaboration, G. von Hippel, R. Sommer, J. Heitger, S. Schaefer and N. Tantalo, 
Ds physics from fine lattices, PoS" LATTICE2008 (2008) 227 [0810.0214]. 

[20] ALPHA Collaboration, M. Delia Morte, R. Hoffmann and R. Sommer, Non-perturbative 
improvement of the axial current for dynamical Wilson fermions, J. High Energy Phys. 03 
(2005) 029 [hep-lat/0503003]. 

[21] T. Bliattacharya, R. Gupta, W.-J. Lee, S. R. Sharpe and J. M. S. Wu, Improved bilinears in 
lattice QCD with non- degenerate quarks, Phys. Rev. D73 (2006) 034504 [hep-lat/0511014]. 

[22] T. Bhattacharya, R. Gupta, W.-J. Lee and S. R. Sharpe, 0(a) improved renormalization 
constants, Phys. Rev. D63 (2001) 074505 [hep-lat/0009038]. 

[23] R. Sommer. Non-perturbative QCD: renormalization, 0(a)-improvement and matching to 
Heavy Quark Effective Theory, Lectures given at ILFTN Workshop on "Perspectives in 
Lattice QCD", Nara, Japan, 31 October -11 November 2005 [hep-lat/0611020]. 

[24] M. Liischcr, S. Sint, R. Sommer and P. Weisz, Chiral symmetry and 0(a) improvement in 
lattice QCD, Nucl. Phys. B478 (1996) 365 [liep-lat/9605038]. 

[25] ALPHA Collaboration, M. Liischer, S. Sint, R. Sommer, P. Wcisz and U. Wolff, 
Non-perturbative 0(a) improvement of lattice QCD, Nucl. Phys. B491 (1997) 323 
[hep-lat/9609035]. 

[26] ALPHA Collaboration, M. Delia Morte, R. Sommer and S. Takcda, On cutoff effects in 
lattice QCD from short to long distances, Phys. Lett. B672 (2009) 407 [0807.1120]. 

[27] ALPHA Collaboration, M. Delia Morte, P. Fritzsch, J. Heitger, H. B. Meyer, H. Simma and 
R. Sommer, Towards a non-perturbative matching of HQET and QCD with dynamical light 
quarks, Po^ LAT2007 (2007) 246 [0710.1188]. 

[28] ALPHA Collaboration, M. Delia Morte, P. Fritzsch, J. Heitger and R. Sommer, 

Non-perturbative quark mass dependence in the heavy-light sector of two-flavour QCD, PoS 
LATTICE2008 (2008) 226 [0810.3166]. 

[29] L. Del Debbio, L. Giusti, M. Liischer, R. Petronzio and N. Tantalo, QCD with light Wilson 
quarks on fine lattices (I): first experiences and physics results, J. High Energy Phys. 02 
(2007) 056 [hep-lat/0610059]. 

[30] ALPHA Collaboration, M. Delia Morte, P. Fritzsch, B. Leder, H. B. Meyer, H. Simma, 
R. Sommer. S. Takeda, O. Witzel and U. Wolff, Preparing for N{ — 2 simulations at small 
lattice spacings, Fo5 LAT2007 (2007) 255 [0710.1263]. 

[31] M. Hasenbusch, Speeding up the Hybrid- Monte- Carlo algorithm for dynamical fermions, 
Phys. Lett. B519 (2001) 177 [hep-lat/0107019]. 

[32] M. Hasenbusch and K. Jansen, Speeding up lattice QCD simulations with clover-improved 
Wilson fermions, Nucl. Phys. B659 (2003) 299 [hep-lat/0211042]. 

[33] ALPHA Collaboration, M. Delia Morte, F. Kncchth, J. Rolf, R. Sommer, I. Wetzorke and 
U. Wolff, Simulating the Schrddinger functional with two pseudo-fermions, Comput. Phys. 
Commun. 156 (2003) 62 [hep-lat/0307008]. 

[34] ALPHA Collaboration, H. B. Meyer, H. Simma, R. Sommer, M. Delia Morte, O. Witzel and 
U. Wolff, Exploring the HMC trajectory-length dependence of autocorrelation times in lattice 
QCD, Comput. Phys. Commun. 176 (2007) 91 [hep-lat/0606004]. 



-29- 



[35] ALPHA Collaboration, U. Wolff, Monte Carlo errors with less errors, Comput. Phys. 
Commun. 156 (2004) 143 [hep-lat/0306017]. 

[36] ALPHA Collaboration, S. Sint and P. Wcisz, Further results on 0(a) improved lattice QCD 
to one-loop order of perturbation theory, Nucl. Phys. B502 (1997) 251 [hep-lat/9704001]. 

[37] ALPHA Collaboration, M. Delia Morte, R. Hoffmann, F. Kneclitli and U. Wolff, Impact of 
large cutoff-effects on algorithms for improved Wilson fermions, Comput. Phys. Commun. 
165 (2005) 49 [hep-lat/0405017]. 

[38] L. D. Dcbbio, L. Giusti, M. Liischer, R. Petronzio and N. Tantalo, Stability of lattice QCD 
simulations and the thermodynamic limit, J. High Energy Phys. 02 (2006) Oil 
[hep-lat/0512021]. 

[39] ALPHA Collaboration, M. Delia Morte, R. Frezzotti, J. Heitger, J. Rolf, R. Sonimcr and 
U. Wolff, Computation of the strong coupling in QCD with two dynamical flavours, Nucl. 
Phys. B713 (2005) 378 [hep-lat/0411025]. 

[40] ALPHA Collaboration, M. Delia Morte, R. Hoffmann, F. Knechtli, J. Rolf, R. Sommer, 
I. Wetzorke and U. Wolff. Non-perturbative quark mass renormalization in two-flavor QCD, 
Nucl. Phys. B729 (2005) 117 [hep-lat/0507035]. 

[41] ALPHA Collaboration, M. Delia Morte, N. Garron, M. Papinutto and R. Sommer, Heavy 
quark effective theory computation of the mass of the bottom quark, J. High Energy Phys. 01 
(2007) 007 [hep-ph/0609294]. 

[42] ALPHA Collaboration, B. Blossier, M. Delia Morte, N. Garron and R. Sommer, HQET at 
order 1/m: I. Non-perturbative parameters in the quenched approximation, 1001.4783. 

[43] ALPHA Collaboration, J. Heitger. A. Jiittner, R. Sommer and J. Wennekers, 

Non-perturbative tests of heavy quark effective theory, J. High Energy Phys. 11 (2004) 048 
[hep-ph/0407227]. 

[44] P. Fritzsch, B-meson properties from non-perturbative matching of HQET to finite-volume 
two-flavour QCD, Ph.D. Thesis (2009). Universitat Miinster. 



-30- 



